1. Open the repo at https://github.com/dlab-berkeley/Geospatial-Fundamentals-in-R-with-sf

   • Clone the repo or download and unzip the zip file
   • Take note of the folder in which the files are located

2. Start RStudio and open a new script, or ./docs/01-core_concepts_and_plotting.Rmd

3. Set your working directory to the folder you unzipped

4. Install the required libraries in RStudio, if you do not have them already

install.packages(c("ggplot2", "sf", "tmap", "raster"),
  dependencies=TRUE)
}

1. Open the slides, ./docs/01-core_concepts_and_plotting.Rmd, in your browser (or click the “Part 1 Slides” link the repo).

Who are you?

Why are you here?

You have options:

1. Open the slides in your browser (click on the link in the repo, or launch from your local copy of the HTML file).

2. Open the Rmd script (or a blank script) in RStudio.

3. Both of the above.

4. Just chill and watch.

Slides ./docs/01-core-concepts-and-plotting.html

RMarkdown Code ./docs/01-core-concepts-and-plotting.Rmd

Make sure you can cut and paste into RStudio

Intro to working with geospatial data in R

• Geospatial data files and formats
• Loading geospatial data in R
• R packages for working with geospatial data
• Coordinate reference systems & map projections
• Mapping geospatial data

Key goal is familiarity, competency takes lots of practice

are data with locations on or near the surface of the Earth.

represent location more specifically with coordinates

46.130479, -117.134167

Coordinates only make sense when associated with a CRS!

Geographic Coordinates: Latitude and Longitude

Define:

• the shape of the Earth

• the origin (0,0 point)

• the relationship between the system and the real world

• the units

Because of variations in these, there are many geographic CRSs!

A CRS can be projected, or can be unprojected (a.k.a. geographic).

A map projection is a mathematial transformation from a curved to a flat surface.

All introduce distortion, in:

• shape
• area
• distance
• direction
• some combination of the above

The larger the area, the greater the distortion!

There is no one map projection best for all purposes. Selection depends on location, extent and purpose.

But some projections are definitely worse than others.

For unprojected data, the World Geodetic System of 1984, or ‘WGS84’, is the most widely used coordinate reference system.

It is the default CRS for most GIS software.

Almost all longitude and latitude assumed to be WGS84 unless otherwise specified.

Historical data much trickier

The ‘connect the dots’ model.

Points, lines and Polygons

We focus on vector data in Parts I and II of this workshop.

The ‘color the squares’ model.

A regular grid of cells (or pixels)

We focus on raster data in Part 3 of this workshop.

• ability to import, create, store, edit, visualize and analyze geospatial data

• representation of geometric data objects, referenced to the surface of the earth via CRSs

• methods to operate on those representations

We call software for working with geospatial data GIS.

This stands for Geographic Information System.

This term is commonly associated with desktop software applications.

Desktop GIS - ArcGIS, QGIS

Spatial Databases - PostgreSQL/PostGIS

Web-based GIS - ArcGIS Online, CARTO

Software geospatial data support - Tableau

Programming languages with geospatial data support

• R, Python, Javascript

You can:

• dynamically determine spatial metrics (e.g. area, length, distance, direction)

• operate on data based on spatial relationships (e.g. intersects, inside, contains)

• link data by location (e.g. linking census data and crime data)

You already use R

Reproducibility

Free & Open Source

Strong support for geospatial data and analysis

Cutting-edge

Make sure you have set your working directory to the location of the workshop files.

Make sure you have installed the packages we are going to use.

#
# Geospatial Data in R Workshop
#

# Make sure needed packages are installed
our_packages<- c("ggplot2", "sf", "tmap", "raster")

for (i in our_packages){
  if ( i %in% rownames(installed.packages()) == FALSE) {
    print(paste(i, "needs to be installed!"))
    install.packages(i)
    
  } else {
    print(paste0(i, " [", packageVersion(i), "] is already installed!"))
    
  } 
}

# Set working directory to folder with workshop files
#setwd("~/Documents/Dlab/workshops/2019/Geospatial-Fundamentals-in-R-with-sf")

There are many approaches to and packages for working with geospatial data in R.

One approach is to keep it simple and store geospatial data in a data frame.

This approach is most common when

• the data are point data in CSV files and

• you want to map rather than spatially transform or analyze the data

sf_properties_25ksample.csv

San Francisco Open Data Portal https://data.sfgov.org

SF Property Tax Rolls

This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.

We are using a subset of these data as a proxy for home values.

SFhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
SFhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

Make sure your working directory is set to the folder where you downloaded the workshop files!

SFhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
SFhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

##   YearBuilt totvalue AreaSquareFeet       Neighborhood NumBedrooms
## 1      1923  1031391           2250 West of Twin Peaks           0
## 2      1909   775300           5830   Presidio Heights           3
## 3      1915  1116772           1350     Inner Richmond           2
## 4      1947   860130           1162    Sunset/Parkside           0
## 5      1981   322749           1463     Outer Richmond           3

class(SFhomes)            # what is the data object type?
dim(SFhomes)              # how many rows and columns
str(SFhomes)              # display the structure of the object
head(SFhomes)             # take a look at the first 10 records
summary(SFhomes)          # explore the range of values
summary(SFhomes$totvalue) # explore the range of values for one column
hist(SFhomes$totvalue)    # histogram for the totvalue column

• What columns contain the geographic data?
• Are these data vector or raster data?
• What type of geometry do the data contain?
  • Points, lines, polygons, grid cells?
• What is the CRS of these data?

We can use the R base plot function to create a simple map

plot(SFhomes$lon, SFhomes$lat) # using base plot function

We can use the R base plot function to create a simple map

plot(SFhomes$lon, SFhomes$lat) # using base plot function

plot(SFhomes$lon, SFhomes$lat) # using base plot function

Let’s add another geospatial data layer to our map.

Map overlay is the fundamental technique of visual geographical analysis.

Use the read.csv function to read in a file of Bart Station locations.

What is the name of the column with the longitude values? latitude?

bart <- read.csv("./data/bart.csv")

# take a look
head (bart)

##           X        Y             STATION OPERATOR DIST  CO
## 1 -122.2833 37.87406      NORTH BERKELEY     BART    4 ALA
## 2 -122.2682 37.86969   DOWNTOWN BERKELEY     BART    4 ALA
## 3 -122.2701 37.85321               ASHBY     BART    4 ALA
## 4 -122.2518 37.84451           ROCKRIDGE     BART    4 ALA
## 5 -122.2671 37.82871           MACARTHUR     BART    4 ALA
## 6 -122.2684 37.80808 19TH STREET/OAKLAND     BART    4 ALA

For the maps from here on out, to deal with a smaller example dataset, we’re going to also subset our data for only those rows that pertain to year 2015.

SFhomes15 <- subset(SFhomes, as.numeric(SalesYear) == 2015)

We can use R’s points function to add points to the current plot.

plot(SFhomes$lon, SFhomes$lat)
points(bart$X, bart$Y, col='red')

Now we’ll add SF landmarks, saved at data/landmarks.csv

landmarks <- read.csv("./data/landmarks.csv")
head(landmarks)

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

Let’s create a map of the SF homes, BART Stations and Landmarks all together.

plot(SFhomes$lon, SFhomes$lat)
points(bart$X, bart$Y, col='red')
points(landmarks$X, landmarks$Y, col='green')

Where’s my third layer? What happened?…

##           X        Y             STATION OPERATOR DIST  CO
## 1 -122.2833 37.87406      NORTH BERKELEY     BART    4 ALA
## 2 -122.2682 37.86969   DOWNTOWN BERKELEY     BART    4 ALA
## 3 -122.2701 37.85321               ASHBY     BART    4 ALA
## 4 -122.2518 37.84451           ROCKRIDGE     BART    4 ALA
## 5 -122.2671 37.82871           MACARTHUR     BART    4 ALA
## 6 -122.2684 37.80808 19TH STREET/OAKLAND     BART    4 ALA

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

Our landmarks CRS seems to be totally different!

How do we reconcile its CRS with the CRS of our other data?…

We cannot yet define and transform CRS.

We cannot yet read and plot other non-CSV vector file formats (such as the very common ESRI Shapefile).

We do not yet know how to answer spatial analysis questions, like:

• What properties are in the Noe Valley neighborhood?

• What is the average property value in each SF neighborhood?

• What is the area of each SF Neighborhood, and the property density?

• What properties are within walking distance (.25 miles) of the Mission neighborhood?

For this, we will need other R packages, with geospatial data structures and functions.

This will allow us to turn R into a fully functional GIS!

Coming up next…

sf stands for ‘simple features’ (a standard developed by the Open Geospatial Consortium).

The standard allows for storing various geometry types in a hierarchical data model.

A ‘feature’ is a representation of a thing in the real world (e.g. a building, a city…).

In other words, each feature consists of both:

• a geometric representation of an object

• some associated information (building: height, name, etc…, city: population, area, etc…).

A relatively new package

supersedes the package sp and its ecosystem (mainly rgeos and rgdal).

A one-stop shop for core geospatial data objects and operations.

sp still frequently used, but its spatial objects are ‘clunkier’.

We’ll get brief exposure to sp in Part 3.

Here are the most common simple features geometries, which are used to represent vector data in sf.

Point Linestring Polygon Multipoint MultiLinestring MultiPolygon Geometry Collection

(From the Geocomputation in R textbook, Chapter 2, Figure 2.2)

sf offers numerous specific benefits, including:

• fast IO (Input and Output)
• enhanced plotting
• integration with R data structures (uses data.frames)
• integration with tidyverse packages and %>% piping syntax
• consistent function names (all start with st_, for ‘spatial type’)
• increasing support by other spatial packages (e.g. tmap, ggplot)
• the simple features data model matches other spatial software (e.g. QGIS, PostGIS, Python’s Geopandas)

First, of course, we’ll need to load the package:

library(sf)

We can then read a spatial dataset into an sf object using the st_read function.

Let’s start with a shapefile of San Francisco census tracts.

First take a look at the file…

dir("data", pattern = "sftracts.")

## [1] "sftracts_wpop.dbf" "sftracts_wpop.prj" "sftracts_wpop.shp"
## [4] "sftracts_wpop.shx" "sftracts.dbf"      "sftracts.prj"     
## [7] "sftracts.shp"      "sftracts.shx"

Now let’s use the sf function st_read to load the file in R.

tracts = st_read(dsn = './data', layer = 'sftracts')

## Reading layer `sftracts' from data source `/home/drew/Desktop/stuff/berk/dlab/Geospatial-Fundamentals-in-R-with-sf/data' using driver `ESRI Shapefile'
## Simple feature collection with 195 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13638250 ymin: 4538276 xmax: -13617440 ymax: 4560139
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs

Then, as always, we can explore the basic aspects of the object returned, using base R functions:

#the object displays a compact summary, when its name is called
tracts

## Simple feature collection with 195 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13638250 ymin: 4538276 xmax: -13617440 ymax: 4560139
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
## First 10 features:
##    STATEFP COUNTYFP TRACTCE             AFFGEOID       GEOID   NAME LSAD
## 1       06      075  035202 1400000US06075035202 06075035202 352.02   CT
## 2       06      075  042601 1400000US06075042601 06075042601 426.01   CT
## 3       06      075  047801 1400000US06075047801 06075047801 478.01   CT
## 4       06      075  060502 1400000US06075060502 06075060502 605.02   CT
## 5       06      075  980200 1400000US06075980200 06075980200   9802   CT
## 6       06      075  018000 1400000US06075018000 06075018000    180   CT
## 7       06      075  012501 1400000US06075012501 06075012501 125.01   CT
## 8       06      075  015200 1400000US06075015200 06075015200    152   CT
## 9       06      075  026403 1400000US06075026403 06075026403 264.03   CT
## 10      06      075  032901 1400000US06075032901 06075032901 329.01   CT
##      ALAND                       geometry
## 1   372758 MULTIPOLYGON (((-13637736 4...
## 2   293133 MULTIPOLYGON (((-13634865 4...
## 3   476606 MULTIPOLYGON (((-13636354 4...
## 4   276334 MULTIPOLYGON (((-13628231 4...
## 5  1022568 MULTIPOLYGON (((-13637838 4...
## 6   937021 MULTIPOLYGON (((-13626895 4...
## 7   135257 MULTIPOLYGON (((-13627070 4...
## 8   279195 MULTIPOLYGON (((-13629456 4...
## 9   379382 MULTIPOLYGON (((-13626769 4...
## 10  671900 MULTIPOLYGON (((-13636086 4...

What sort of object is this?

class(tracts)

class(tracts)

## [1] "sf"         "data.frame"

The object is of both the ‘sf’ and ‘data.frame’ classes.

It has a number of columns (i.e. attributes, fields), including a geometry column

str(tracts)

## Classes 'sf' and 'data.frame':   195 obs. of  9 variables:
##  $ STATEFP : Factor w/ 1 level "06": 1 1 1 1 1 1 1 1 1 1 ...
##  $ COUNTYFP: Factor w/ 1 level "075": 1 1 1 1 1 1 1 1 1 1 ...
##  $ TRACTCE : Factor w/ 195 levels "010100","010200",..: 164 169 178 184 191 70 26 42 129 155 ...
##  $ AFFGEOID: Factor w/ 195 levels "1400000US06075010100",..: 164 169 178 184 191 70 26 42 129 155 ...
##  $ GEOID   : Factor w/ 195 levels "06075010100",..: 164 169 178 184 191 70 26 42 129 155 ...
##  $ NAME    : Factor w/ 195 levels "101","102","103",..: 164 169 178 184 191 70 26 42 129 155 ...
##  $ LSAD    : Factor w/ 1 level "CT": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ALAND   : num  372758 293133 476606 276334 1022568 ...
##  $ geometry:sfc_MULTIPOLYGON of length 195; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:9, 1:2] -13637736 -13636969 -13636954 -13636940 -13636925 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr  "STATEFP" "COUNTYFP" "TRACTCE" "AFFGEOID" ...

We can use basic data.frame functions on it

nrow(tracts)

## [1] 195

colnames(tracts)

## [1] "STATEFP"  "COUNTYFP" "TRACTCE"  "AFFGEOID" "GEOID"    "NAME"    
## [7] "LSAD"     "ALAND"    "geometry"

head(tracts, 4)

## Simple feature collection with 4 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13637740 ymin: 4538290 xmax: -13627250 ymax: 4549243
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
##   STATEFP COUNTYFP TRACTCE             AFFGEOID       GEOID   NAME LSAD
## 1      06      075  035202 1400000US06075035202 06075035202 352.02   CT
## 2      06      075  042601 1400000US06075042601 06075042601 426.01   CT
## 3      06      075  047801 1400000US06075047801 06075047801 478.01   CT
## 4      06      075  060502 1400000US06075060502 06075060502 605.02   CT
##    ALAND                       geometry
## 1 372758 MULTIPOLYGON (((-13637736 4...
## 2 293133 MULTIPOLYGON (((-13634865 4...
## 3 476606 MULTIPOLYGON (((-13636354 4...
## 4 276334 MULTIPOLYGON (((-13628231 4...

We can plot an sf object using its plot method.

In other words, when we just call R’s base plot function on an sf object, R will recognize that it’s an sf object and thus plot it accordingly.

plot(tracts)

plot(tracts)

We get an array of plots, one for each variable (up to the first 9).

So then we should be able to plot a single variable by just plotting a subsetted sf dataframe.

Plot just the ‘NAME’ column’s data.

(Note: This will be an example of what we call a ‘choropleth’ map.)

#read in a shapefile of SF census tracts
plot(tracts['NAME'])

Some of you may have gotten this plot instead:

plot(tracts$NAME)

What went wrong?

class(tracts['NAME'])

## [1] "sf"         "data.frame"

class(tracts[, 'NAME'])

## [1] "sf"         "data.frame"

class(subset(tracts, select='NAME'))

## [1] "sf"         "data.frame"

class(tracts$NAME)

## [1] "factor"

When we use bracket syntax or the subset function, sf objects return new, subsetted sf objects.

But when we use the ‘$’ notation, we just get a vector of the column’s values!

As always, we need to be careful and check what kinds of objects we’re working with!

head(tracts$NAME)

## [1] 352.02 426.01 478.01 605.02 9802   180   
## 195 Levels: 101 102 103 104 105 106 107 108 109 110 111 112 113 117 ... 9809

head(tracts['NAME'])

## Simple feature collection with 6 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13637840 ymin: 4538290 xmax: -13624700 ymax: 4549564
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
##     NAME                       geometry
## 1 352.02 MULTIPOLYGON (((-13637736 4...
## 2 426.01 MULTIPOLYGON (((-13634865 4...
## 3 478.01 MULTIPOLYGON (((-13636354 4...
## 4 605.02 MULTIPOLYGON (((-13628231 4...
## 5   9802 MULTIPOLYGON (((-13637838 4...
## 6    180 MULTIPOLYGON (((-13626895 4...

The nice thing about sf objects is that they are just data.frames!

The geometry data is just stored in its own special column, usually named geom or geometry.

For the most part, we will not want to manually manipulate the data in the geometry column.

But when we’re just getting started, it can be revealing to tinker a bit.

As we saw earlier, the geometry data is stored in its own column. Let’s take a closer look at that column:

tracts$geometry

## Geometry set for 195 features 
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13638250 ymin: 4538276 xmax: -13617440 ymax: 4560139
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
## First 5 geometries:

## MULTIPOLYGON (((-13637736 4546153, -13636969 45...

## MULTIPOLYGON (((-13634865 4549210, -13634268 45...

## MULTIPOLYGON (((-13636354 4548053, -13635877 45...

## MULTIPOLYGON (((-13628231 4538555, -13628045 45...

## MULTIPOLYGON (((-13637838 4548684, -13637471 45...

We see that our CRS and some other metadata are connected with the data in this column.

This includes:

• the geometry type and its dimensionality
• the bounding box,
• the CRS arguments
• the CRS’ EPSG code (which we’ll learn about in a bit)

tracts$geometry

## Geometry set for 195 features 
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13638250 ymin: 4538276 xmax: -13617440 ymax: 4560139
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
## First 5 geometries:

## MULTIPOLYGON (((-13637736 4546153, -13636969 45...

## MULTIPOLYGON (((-13634865 4549210, -13634268 45...

## MULTIPOLYGON (((-13636354 4548053, -13635877 45...

## MULTIPOLYGON (((-13628231 4538555, -13628045 45...

## MULTIPOLYGON (((-13637838 4548684, -13637471 45...

All the geometries are always in the same CRS.

We can check that CRS using the st_crs function:

st_crs(tracts)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs"

We can also see that this column is some sort of special object.

Unlike with the ‘NAME’ column, when we subset the ‘geometry’ column with ‘$’ we don’t just get a vector of values.

This is because we can’t just get a vector, because the values are complex objects (geometries), not single values!

So instead, we get a “list column”.

So what sort of object is the geometry column?

class(tracts$geometry)

## [1] "sfc_MULTIPOLYGON" "sfc"

class (tracts$geometry[[1]])

## [1] "XY"           "MULTIPOLYGON" "sfg"

The column is an sfc object (‘simple features collection’).

It is composed of sfg objects (‘simple features geometry’).

For our purposes, we just need to think of this as a list of geometries!

We will see each geometry (sfg) in our column represented as well-known text (WKT).

Here are some examples:

POINT (0 3)

LINESTRING (0 3, 1 4, 2 3)

POLYGON ((5 2, 1 7, 1 1))

These are all the dots in the “connect the dots” model!

That is the full anatomy of an sf object!

• sf objects are data.frame objects with special ‘geom’ or ‘geometry’ columns.
• The ‘geometry’ column is a simple features collection (sfc) object.
• Each value in an sfc is a simple features geometry (sfg) object.

(From the sf docs)

Here’s a more technical summary (from sf’s own Github Repo):

Now if we look at our census tracts again, we can see all the pieces we just reviewed:

tracts

## Simple feature collection with 195 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -13638250 ymin: 4538276 xmax: -13617440 ymax: 4560139
## epsg (SRID):    NA
## proj4string:    +proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs
## First 10 features:
##    STATEFP COUNTYFP TRACTCE             AFFGEOID       GEOID   NAME LSAD
## 1       06      075  035202 1400000US06075035202 06075035202 352.02   CT
## 2       06      075  042601 1400000US06075042601 06075042601 426.01   CT
## 3       06      075  047801 1400000US06075047801 06075047801 478.01   CT
## 4       06      075  060502 1400000US06075060502 06075060502 605.02   CT
## 5       06      075  980200 1400000US06075980200 06075980200   9802   CT
## 6       06      075  018000 1400000US06075018000 06075018000    180   CT
## 7       06      075  012501 1400000US06075012501 06075012501 125.01   CT
## 8       06      075  015200 1400000US06075015200 06075015200    152   CT
## 9       06      075  026403 1400000US06075026403 06075026403 264.03   CT
## 10      06      075  032901 1400000US06075032901 06075032901 329.01   CT
##      ALAND                       geometry
## 1   372758 MULTIPOLYGON (((-13637736 4...
## 2   293133 MULTIPOLYGON (((-13634865 4...
## 3   476606 MULTIPOLYGON (((-13636354 4...
## 4   276334 MULTIPOLYGON (((-13628231 4...
## 5  1022568 MULTIPOLYGON (((-13637838 4...
## 6   937021 MULTIPOLYGON (((-13626895 4...
## 7   135257 MULTIPOLYGON (((-13627070 4...
## 8   279195 MULTIPOLYGON (((-13629456 4...
## 9   379382 MULTIPOLYGON (((-13626769 4...
## 10  671900 MULTIPOLYGON (((-13636086 4...

So what about our homes data? Let’s plot them with the tracts data!

NOTE: With the base plot function, we need to use sf::st_geometry to pull out just the sf object’s geometric features, for plotting.

plot(sf::st_geometry(tracts))
points(SFhomes15$lon, SFhomes15$lat)

What happened? Where are the homes?

Once again, we cannot plot our data together because the CRS’s are not the same!

st_crs(tracts)

## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +R=6378137 +units=m +no_defs"

st_crs(SFhomes15)

## Coordinate Reference System: NA

The #1 reason…

Every spatial object needs a reference to its CRS.

Because the homes are lacking this reference, we’ll need to define it.

Note the important difference between defining and transforming projections:

• Defining just makes a spatial object aware of the CRS of its geometries. (This does not change the data’s coordinate values.)

• Transforming converts the geometries to a new CRS. (This changes the coordinate values.)

So we need our data all in the same projection.

In order to do this, we’re going to need to follow 3 steps:

1. Coerce the homes to an explicitly geospatial object, using sf.

2. Define (or assign, or set) that object’s CRS.

3. Transform one of our two objects to the CRS of the other.

With sf, we can combine steps 1 and 2 in a single step, using the st_as_sf function.

To do this, we need to provide 3 pieces of information to 3 arguments:

x: The object to be converted to an sf object

coords: The columns containing the coordinates.

crs: The CRS of the coordinates contained in those columns.

In order to this, we need to tell R the CRS of our homes data.

What is our CRS?

And how can well tell that to R?

GIS software has database of definitions for thousands of Earth-referenced CRSs.

Those definitions are are used by GIS, to transform data between CRS’s.

We just need to know how to specify our data’s CRS in terms that will allow R to correctly indentify it within that database.

Our data’s CRS is unprojected longitude/latitude on a WGS84 datum.

We can specify it using a common format for referencing CRS’s, the proj4 string:

st_crs("+proj=longlat +ellps=WGS84 +datum=WGS84")

However, the proj4 string is long and complicated, and leaves a lot of room for error.

So the European Petroluem Survey Group devised a simpler system of numerical codes, called EPSG codes.

The longlat, WGS84 EPSG code is 4326.

So we could create our CRS instead using just the EPSG code in a proj4string:

st_crs("+init=epsg:4326")

Which, in sf, can be simplified even further:

st_crs(4326)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Here is some additional info on CRS definitions, for your reference:

Proj4 is the standard library for defining and transforming map projections and coordinate reference systems.

Here is an example file of proj4 CRS definitions

And here are some common codes:

Geographic CRSs

• 4326 Geographic, WGS84 (default for lon/lat)

• 4269 Geographic, NAD83 (USA Fed agencies like Census)

Projected CRSs

• 5070 USA Contiguous Albers Equal Area Conic

• 3310 CA ALbers Equal Area

• 26910 UTM Zone 10, NAD83 (Northern Cal)

• 3857 Web Mercator (web maps)

And here’s a resource for looking up EPSG codes and proj4 strings, should you need to:

http://spatialreference.org/

Now that we know our homes have are in unprojected coordinates with a WGS84 datum (i.e. EPSG code: 4326), we can call st_as_sf as follows:

SFhomes15_sf = st_as_sf(SFhomes15, coords = c('lon', 'lat'), crs = 4326)

We can check the object type by displaying a summary.

SFhomes15_sf

## Simple feature collection with 835 features and 16 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -122.5107 ymin: 37.70854 xmax: -122.3696 ymax: 37.80612
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##     FiscalYear  SalesDate                               Address YearBuilt
## 24        2015 2015-08-21  0000 2760 19TH                AV0015      1979
## 35        2015 2015-08-13  0000 0560AMISSOURI            ST0000      2003
## 76        2015 2015-12-29  0000 0718 LONG BRIDGE         ST1202      2016
## 79        2015 2015-08-26  0000 0103 SUMMIT              WY0000        NA
## 85        2015 2015-07-06  0000 0899 VALENCIA            ST0202      2015
## 90        2015 2015-06-12  0000 1333 JONES               ST0808      1961
## 92        2015 2015-10-08  0000 0167 TOPAZ               WY0000      1973
## 98        2015 2015-04-14  0000 1904 BAKER               ST0000      1900
## 223       2015 2015-07-23  0000 0301 SUMMIT              WY0000        NA
## 229       2015 2015-09-18 0201 0201 FOLSOM              ST00A3K      2015
##     NumBedrooms NumBathrooms NumRooms NumStories NumUnits AreaSquareFeet
## 24            2            2        5          0        1           1595
## 35            2            2        5          1        1           1191
## 76            2            2        5          1        1           1346
## 79            3            3        6          0        1           2133
## 85            3            3        7          1        1           1266
## 90            3            3        7         14        1           1840
## 92            3            3        7          2        1           2075
## 98            2            2        5          1        1           1256
## 223           0            0        0          0        0              0
## 229           2            2        0          1        1           1520
##     ImprovementValue LandValue                   Neighborhood SupeDistrict
## 24            432500    432500             West of Twin Peaks            7
## 35            701280    701280                   Potrero Hill           10
## 76           1512093    748900                    Mission Bay            6
## 79            582108    873161                      Lakeshore            7
## 85            850000    850000                        Mission            9
## 90           1154846   1154846                       Nob Hill            3
## 92            484544   1130604                      Glen Park            8
## 98            810169   1890395                Pacific Heights            2
## 223           486872   1136036                      Lakeshore            7
## 229           962500    962500 Financial District/South Beach            6
##     totvalue SalesYear                   geometry
## 24    865000      2015 POINT (-122.4741 37.73601)
## 35   1402560      2015  POINT (-122.3965 37.7592)
## 76   2260993      2015 POINT (-122.3942 37.77181)
## 79   1455269      2015 POINT (-122.4754 37.71518)
## 85   1700000      2015  POINT (-122.421 37.75876)
## 90   2309692      2015 POINT (-122.4149 37.79362)
## 92   1615148      2015   POINT (-122.437 37.7402)
## 98   2700564      2015 POINT (-122.4437 37.78881)
## 223  1622908      2015 POINT (-122.4754 37.71492)
## 229  1925000      2015 POINT (-122.3919 37.78868)

We can also specifically check the CRS with st_crs:

st_crs(SFhomes15_sf)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Looks good!

That takes care of both steps 1 and 2!

Now for step 3.

We need to take one data set from its current CRS to the other data set’s CRS, to reconcile the two.

This is called transforming, or reprojecting.

We want all data in the same CRS.

Thus, we can either transform SFhomes15_sf to the tracts’s CRS, or vice versa.

Which one is best???

There’s no universally right answer.

Usually depends on what operations we plan to run downstream.

For now, we’ll reproject the tracts.

That way, everything will be in an unprojected CRS, and will plot in decimal degrees (rather than meters, for example).

Note: If we were working with sp objects, we would use the rgdal package to carry out CRS transformations.

However, sf encapsulates rgdal’s projection-related capabilities.

Using sf, we can reproject, or transform an sf object using the st_transform function.

The only arguments we need are:

• x: The object to be transformed

• crs: The target CRS to transform x to (which can be provided as just an integer of the EPSG code)

Syntax…

new_sf_object = st_transform(original_sf_object, crs=target_crs_epsg_code)

Transform the tracts object to the same CRS as the SFhomes15_sf object.

Name the new object tracts_lonlat (to indicate that its CRS is unprojected longitude and latitude).

We can transform the tracts to the CRS of SFhomes15_sf using:

tracts_lonlat = st_transform(tracts, crs = 4326)

However, note that we could be more explicit!

This approach will make certain that things match up:

tracts_lonlat = st_transform(tracts, crs = st_crs(SFhomes15_sf))

Are they both defined? Are they the same?

st_crs(SFhomes15_sf)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

st_crs(tracts_lonlat)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

st_crs(SFhomes15_sf) == st_crs(tracts_lonlat)

## [1] TRUE

And now we can try again to plot our data together.

(Note: Now that both objects are sf objects, we can actually just call the plot method for each, as long as we pass it the add=TRUE argument.)

plot(st_geometry(tracts_lonlat))
plot(st_geometry(SFhomes15_sf), add=TRUE, col='blue')

Success!

Work through the following steps:

1.) Go back to the landmarks data that we stranded. Coerce it from a dataframe to an sf object, setting its CRS to EPSG: 3857, which is the code for Web Mercator.

2.) Read in the ‘SFboundary’ and ‘SFhighways’ shapefiles, from the ‘./data’ subdirectory.

3.) Reproject either of those layers to the CRS of SFhomes15_sf, as needed.

4.) Plot tracts, boundary, highways, homes, and landmarks together. Make the boundary’s border purple (check the docs!).

5.) Make the highways and the landmarks red, and the homes blue.

1.) Go back to the landmarks data that we stranded. Coerce it from a dataframe to an sf object, setting its CRS to EPSG: 3857, which is the code for Web Mercator.

landmarks_sf = st_as_sf(landmarks, coords = c('X', 'Y'), crs = 3857)

2.) Read in the ‘SFboundary’ shapefile, from the ‘./data’ subdirectory.

SFboundary = st_read('./data', 'SFboundary')

## Reading layer `SFboundary' from data source `/home/drew/Desktop/stuff/berk/dlab/Geospatial-Fundamentals-in-R-with-sf/data' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -122.5151 ymin: 37.70825 xmax: -122.357 ymax: 37.81176
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

SFhighways = st_read('./data', 'SFhighways')

## Reading layer `SFhighways' from data source `/home/drew/Desktop/stuff/berk/dlab/Geospatial-Fundamentals-in-R-with-sf/data' using driver `ESRI Shapefile'
## Simple feature collection with 246 features and 5 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 543789.3 ymin: 4173548 xmax: 551295.9 ymax: 4183929
## epsg (SRID):    26910
## proj4string:    +proj=utm +zone=10 +datum=NAD83 +units=m +no_defs

3.) Reproject either of those layers to the CRS of SFhomes15_sf, as needed.

SFboundary doesn’t need to be transformed.

#check the CRS of SFboundary
st_crs(SFboundary) == st_crs(SFhomes15_sf)

## [1] TRUE

3.) Reproject either of those layers to the CRS of SFhomes15_sf, as needed.

SFhighways needs to be transformed.

#check th CRS of SFhighways
st_crs(SFhighways) == st_crs(SFhomes15_sf)

## [1] FALSE

#it needs to be transformed
SFhighways_lonlat = st_transform(SFhighways, st_crs(SFhomes15_sf))

3.) Reproject either of those layers to the CRS of SFhomes15_sf, as needed.

We know landmarks needs to be reprojected, because we just assinged it EPSG 3857.

landmarks_lonlat = st_transform(landmarks_sf, st_crs(SFhomes15_sf))

4.) Plot tracts, boundary, highways, homes, and landmarks together. Make the boundary’s border purple (check the docs!).

5.) Make the highways and the landmarks red, and the homes blue.

plot(st_geometry(SFboundary), col='white', border='purple')
plot(st_geometry(tracts_lonlat), alpha=0.2, add=T)
plot(st_geometry(SFhighways_lonlat), col='red', add=T)
plot(st_geometry(SFhomes15_sf), col='blue', pch=16, add=T)
plot(st_geometry(landmarks_lonlat), col='green', pch=16, cex=2, add=T)

Excellent! Everything lines up!

How can we manipulate the aspects of our plots?

How can we make our plots prettier?

What other options do we have for plotting sf objects?

We need specific plotting packages!

Most widely used plotting library in R

Not specifically for geospatial data, but can be used to make fabulous maps

Important starting point for the syntax of our go-to mapping package, tmap

Load the library

library(ggplot2)

Here’s how we could plot our homes, colored by their total values:

ggplot() + geom_point(data=SFhomes15, aes(x=lon, y=lat, col=totvalue))

Break down the syntax:

• ggplot(): creates a ggplot object
• ` + geom_point`: adds a point layer to our map
• data=SFhomes15_sf: tells ggplot2 to draw data from the SFhomes15_sf object
• aes: ‘aesthetically maps’ our variables onto aspects of our plot
• x and y: tell ggplot2 which columns to use for the x- and y-axis positions of the points
• col='totvalue': tells ggplot2 to color our points by the values in the ‘totvalue’ column

Data-driven symbology!

ggplot() + geom_point(data=SFhomes15, aes(x=lon, y=lat, col=totvalue))

Using ggplot2 gave us more flexbile symbology than using base R’s plot.

Using tmap will give use more cartographic power than using ggplot2.

We’ll see a quick overview of the package today.

Then we’ll get continued practice with tmap in Parts 2 and 3.

As we will see, tmap offers some crucial additional options:

• integration with spatial objects (i.e. sf and sp objects)
• easy interactive maps

tmap stands for thematic map

It’s a powerful toolkit for creating maps with sf and sp objects, with less code than the alternatives

Syntax inspired by ggplot2, but a bit simpler

First we’l load the library

library(tmap)

## Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.8.0-CAPI-1.13.1
## and GEOS at installation 3.7.2-CAPI-1.11.2differ

Note: You may need to install /load dependencies - ggplot2, RColorbrewer, classInt, leaflet libraries

We can create quick tmaps with the qtm function:

qtm(SFhomes15_sf)

tmap has 2 modes:

• plot: makes static maps

• view: makes interactive maps

Use the command tmap_mode to toggle between them:

• tmap_mode("plot")

• tmap_mode("view")

We can create an interactive version of the previous map:

tmap_mode("view") # set tmap to interactive view mode
qtm(SFhomes15_sf) # Interactive - click on the points

Notice that this is a live, interactive map! Click to see Popups! Click on the Layer Selector!

Conveniently, we can also toggle between the modes using the ttm function.

ttm()

## tmap mode set to plotting

ttm()

## tmap mode set to interactive viewing

ttm()

## tmap mode set to plotting

ttm()

## tmap mode set to interactive viewing

ttm()

## tmap mode set to plotting

To customize our tmaps, we need to use ggplot2-like syntax.

Here’s how we could recreate the plot we just made in ggplot2.

tm_shape(SFhomes15_sf) + tm_dots(col="totvalue", size = 0.25)

Break down the syntax:

• tm_shape(): creates a tmap object
• SFhomes15_sf: tells tmap to draw data from the SFhomes15_sf object
• ` + tm_dots`: adds a dots layer to our map
• col='totvalue': tells tmap to color our dots by their values in the ‘totvalue’ column

tm_shape(SFhomes15_sf) + tm_dots(col="totvalue", size = 0.25)

Like ggplot, we used data values to drive the symbology.

But note three important differences:

1. We didn’t use an aesthetic-mapping function like aes.
2. We didn’t need to specify the x- and y-axis positions.
3. We had to ‘quote’ our column names!

The tm_dots functions added dots to our map.

tm_dots are a type of tm_symbols()

See ?tm_symbols for other types of point symbols.

We can also add and customize polygons using tm_polygons.

tm_shape(tracts) + tm_polygons(col="beige", border.col="red", alpha = 0.4)

And we can add and customize lines using tm_lines.

tm_shape(SFhighways_lonlat) + tm_lines(col="black")

We can overlay multiple sf objects, or layers, by concatenating the tmap geometry-layer commands with plus signs (again, using ggplot-style syntax).

# Map the SF Boundary first
overlay_map = tm_shape(SFboundary) + 
  tm_polygons(col="beige", border.col="black") +

# Overlay the highway lines next
tm_shape(SFhighways_lonlat) + 
  tm_lines(col="black") +
  
# Then add the house points
tm_shape(SFhomes15_sf) + 
  tm_dots(col="totvalue", size=.25) 

overlay_map

What if we then wanted to add our landmarks?

Does this code work?

overlay_map +
    tm_shape(landmarks_lonlat) +
    tm_dots(col = 'skyblue', size = 2)

Yup it does!

What if we redid our map using landmarks_sf instead of landmarks_lonlat?

Does this code work?

overlay_map +
    tm_shape(landmarks_sf) +
    tm_dots(col = 'skyblue', size = 2)

Yup! What does that tell us about tmap?

tmap reprojects on the fly!

• If the CRSs are defined, tmap will use that info.

• If not, tmap will assume WGS84.

• It will dynamically reproject subsequent layers to match the first one.

What’s happening here?

SFhomes15_sf_low2high <- SFhomes15_sf[order(SFhomes15_sf$totvalue, decreasing = FALSE),]
SFhomes15_sf_high2low <- SFhomes15_sf[order(SFhomes15_sf$totvalue, decreasing = TRUE),]

high2low = tm_shape(SFhomes15_sf_high2low) + 
  tm_dots(col='totvalue', size=1)
low2high = tm_shape(SFhomes15_sf_low2high) + 
  tm_dots(col='totvalue', size=1)

tmap_arrange(high2low, low2high)

The order of the data in the data frame changes the map display. This is an important cartographic concern!

We can of course tweak all sorts of details in our maps.

tm_shape(SFboundary) +
  tm_polygons(col="beige", border.col="black") +
  
tm_shape(SFhighways_lonlat) + 
  tm_lines(col="black") +
  
tm_shape(SFhomes15_sf) + 
  tm_dots(col="totvalue", size=.25, 
          title = "San Francisco Property Values (2015)") + 
tm_layout(inner.margins=c(.05, .2, .15, .05)) 
      # bottom, left, top, right

We can also customize the popups!

Notice the changes relative to the previous code.

#set mode back to interactive
tmap_mode('view')

tm_shape(SFboundary) + 
  tm_polygons(col="beige", border.col="black") +
  
tm_shape(SFhighways_lonlat) + 
  tm_lines(col="black") +

tm_shape(SFhomes15_sf) + 
  tm_dots(col="totvalue", size=.05, 
          title = "San Francisco Property Values (2015)",
          popup.vars=c("SalesYear","totvalue","NumBedrooms",
                       "NumBathrooms","AreaSquareFeet")) + 
  
tm_layout(inner.margins=c(.05, .2, .15, .05)) # bottom, left, top, right

And we can even access a wide variety of basemaps.

Here’s the list: http://leaflet-extras.github.io/leaflet-providers/preview/

Let’s try watercolor!

tm_basemap("Stamen.Watercolor") +
tm_shape(SFhomes15_sf) + 
tm_dots(col="totvalue", size=.05, title = "San Francisco Property Values (2015)") +
tm_tiles("Stamen.TonerLabels")

Let’s try watercolor!

We can use the tmap_last function to grab the last map we made.

Here, we save it to a named variable (‘map1’), then display it.

fav_map <- tmap_last()

Then we can display it.

fav_map

And then we can use tmap_save to save it.

Let’s save in a couple formats, then look at them.

Note: We can specify the output format by just using the file extension!

tmap_save(fav_map, "./output/SF_properties.png", height=6) # Static image file with
tmap_save(fav_map, "./output/SF_properties.html") # interactive web map

There are many ways to publish your maps.

You can share you map online by publishing it to RPubs.

• You need to have an RPubs account to make that work.

1. Enter the name of your tmap object (map1) or your tmap code in the console

2. In the Viewer window, click on the Publish icon.

Here is an RPubs Demo.

tmap provides support for so many different types of maps, check out the tmap: get started guide, vignette("tmap-getstarted"), linked in the documentation.

#Check out the tmap: get started! guide
?tmap

You also plenty of options for further information:

• ?tmap

• ?tmap_shape

• ?tmap_element

• ?tm_polygons (tm_fill, tm_borders)

• ?tm_symbols (tm_dots, etc…)

• The Geocomputation with R textbook (Lovelace, Nowosad, and Muenchow, 2019).

Using tmap instead of base plot, recreate the map from our previous challenge, but this time color the homes by their total property values!.

For reference, here’s a summary of the prompt, along with the previous map:

• Plot tracts, boundary, highways, homes, and landmarks together.

• Make the boundary’s border purple.

• Make the highways and the landmarks red

• Color the homes by their total property values.

challenge_map_2 = tm_shape(SFboundary) +
    tm_polygons(border.col = 'purple', alpha = 0) +
tm_shape(tracts_lonlat) +
    tm_polygons(col = 'lightgray', border.col = 'darkgray', alpha = 0.3) + 
tm_shape(SFhighways_lonlat) +
    tm_lines(col = 'red') +
tm_shape(SFhomes15_sf) +
    tm_dots(col = 'totvalue', size = 0.25, palette = '-Blues') +
tm_shape(landmarks_sf) +
    tm_dots(col = 'green', size = 0.5) 

tmap_mode('plot')
challenge_map_2 + tm_layout(inner.margins=c(.05, .2, .15, .05)) 

• Foundations of GIS
• Geospatial Data in R
• Basic plots and overlays
• sf library: spatial objects and methods in R
• Converting data frames to sf
• CRS definitions and transformations
• tmap: advanced mapping

That’s the end of Part I!

See you in Part II, where we’ll cover spatial analysis!